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Chaotic micromixers such as the staggered herringbone mixer developed by Stroock et aL aUow efficient 
mixing of fluids even at low Reynolds number by repeated stretching and folding of the fluid interfaces. 
The ability of the fluid to mix well depends on the rate at which "chaotic advection" occurs in the mixer. 
An optimization of mixer geometries is a non trivial task which is often performed by time consuming and 
expensive trial and error experiments. In this paper an algorithm is presented that applies the concept of 
finite-time Lyapunov exponents to obtain a quantitative measure of the chaotic advection of the flow and 
hence the performance of micromixers. By performing lattice Boltzmann simulations of the flow inside a 
mixer geometry, introducing massless and non-interacting tracer particles and following their trajectories 
the finite time Lyapunov exponents can be calculated. The applicability of the method is demonstrated by a 
comparison of the improved geometrical structure of the staggered herringbone mixer with available literature 
data. 
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I. INTRODUCTION 

Microfluidic devices have found applications in various 
scientific and industrial processes. A typical example is 
their integration as irnportant components of chemical 
and biological sensors.'^' A micromixer is a microfiuidic 
device used for effective mixing of different fluid con- 
stituents. It can be used to efficiently mix for exam- 
ple a variety of bio-reactants such as bacteria cells, large 
DNA molecules, enzymes and proteins in portable in- 
tegrated microsystems with minimum energy consump- 
tion. It is also used in mixing of solutions in chemical 
reactions,'^ sequencing of nucleic acids or drug solution 
dilution. Hence, the design of practical and efficient mi- 
cromixers is a major research topic in microfluidics,'^ 
especially in the development of micro total analysis sys- 
tems. Over the years various methods of efficient mixing 
have been developed and many of those have been suc- 
cessfully applied in industry.^ 

The small length scales of the micromixers has a nega- 
tive impact on mixing as it results in laminar flows inside 
the channels. In this flow regime, mixing is influenced 
mainly by the process of molecular intcr-diffusion.I^ Ex- 
periments with channels with complex surface topology 
including grooved walls have revealed that microscale 
mixing is enhanced by "chaotic advection" , a process 
which was flrst reviewed by Aref in 1984.'^ He describes 
how mixing is still possible even at low Reynolds number 
by repeated stretching and folding of fluid elements.'^ If 
properly applied, this mechanism causes the interfacial 
area between the fluids to increase exponentially, which 
can then lead to an enhanced intermaterial transport, 
hence mixing. A comprehensive mathematical descrip- 



tion of the exponential growth of interfacial surfaces can 
be found in the book by Ottino.^ Mixers were designed in 
the following years, which utilize the principle of "chaotic 
advection" .ES However, it is important to note that the 
term "chaotic" is used strictly in a Lagrangian sense.'^ 

Depending on the working principle, micromixers can 
be categorized into two important types: if energy from 
an external source is used to drive the mixing process, 
then it is termed as "active mixer". These external en- 
ergy sources can be acoustic bubble induced vibrations, 
periodic variation of the flow rate, piezoelectric vibrating 
membranes, valves, etc. The external sources are often 
moving components such as micropumps and they re- 
quire advanced fabrication steps.l^ The second category 
of micromixers is based on restructuring the flow profile 
using static but sophisticated mixer geometries. These 
are termed as "passive mixer" . Passive micromixers have 
the advantage of simple fabrication, easy operation and 
no elements which can generate heat. The absence of 
heating is an important factor for applications to biolog- 
ical studies where temperature is a sensitive parameter. 
The mixing length and mixing time are defined as the dis- 
tance and time span the fluid constituents have to flow in- 
side the mixer in order to obtain a homogeneous mixture. 
An effective micromixer should reduce the mixing length 
and mixing time substantially in order to achieve rapid 
mixing. A common practice to design passive micromix- 
ers is to create alternating thin fluid lamellae. These 
result in an interfacial area that increases with the num- 
ber of lamellae rendering the diffusion process more ef- 
fective and hence allowing fast mixing. 131 There are many 
examples of bi-lamellatiorP^ESI ^^d multi-lamellatiorP in 
the literature on micromixers. However, the drawback 
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of such devices is that the number of lamellae is gener- 
ally limited due to the negative impact on the applied 
pressure drop caused by the required microstructures in- 
side the channel. Other examples of passive micromix- 
ers include the twisted pipe mixer,^ the superfocus mi- 
cromixer, where several jets are made to collide at the 
focal point of the jets and the three dimensional serpen- 
tine modelP 

The recently developed so-called "chaotic micromixer" 
has gained substantial interest in the literature since it 
overcomes some of the drawbacks of conventional mix- 
ers based on muti-lamellation. Such a device consists 
of microstructured objects such as "herringbones" (see 
Fig. [T]). The staggered herringbone mixer (SHM) was 
introduced as one of the first experimental implementa- 
tions of chaotic micromixers in 2002 by Stroock et 
The half cycles of the SHM consist of grooves with two 
arms which are asymmetric and unequal in length. These 
arms are inclined at an angle of 45° to the wall and 90° 
against each other, while the pattern interchanges every 
half cycle of the herringbone. The peculiar arrangement 
of the herringbone structure enhances the mixing process 
by "chaotic advection" where the interfacial area between 
the fluids grows exponentially in time - another impor- 
tant advantage over mixers using the concept of multi- 
lamellation. The SHM was used by several authors for 
studies related to mixing and analysis of mixing quality. 
To characterize the mixing quality of the SHM, Aubin et 
al.ESI implemented particle tracking techniques using the 
variance of tracer dispersion. Li and Chen report on an 
optimization of the SHM using the standard deviation 
of particle concentrations. As in the current paper, they 
used the lattice Boltzmann (LB) method to model the 
flow field.'i^! Further, the SHM was used by Kirtland et 
al^ to study the mass transfer to reactive boundaries 
from three-dimensional micro-channels. 

For the development of micromixers it is important to 
have reliable tools at hand to quantify their performance. 
Efficiency and mixing quality have been studied by var- 
ious methods in the past. These include the analysis 
of the probability density function of the flow proflles, 
studying the stretching of the flow field, the Poincare 
section analysis, or the intensity of segregation as intro- 
duced by Danckwerts in 1952.'2ll In this paper an alter- 
native numerical procedure is presented which is tailored 
for the optimization of chaotic micromixers. It is based 
on LB simulations to describe the flow inside complex 
mixer geometries together with a measurement of flnite 
time Lyapunov exponents (FTLE) as obtained from tra- 
jectories of massless tracer particles immersed in the flow. 
The LB method can easily handle flows in complex ge- 
ometries, which makes this method convenient for flows 
in advanced microstructures such as the micromixers the 
current paper focuses on. The Lyapunov exponent pro- 
vides a quantitative measure of long term average growth 
rates of small initial flow perturbations and thus all ows a 
quantification of the efficiency of chaotic advection .1221221 
Since the systems of interest are finite and simulations are 



limited to a flnite time span, the proposed method uti- 
lizes Wolf's method to calculate the FTLE.i^ While LB 
simulations of flows in micromixers and FTLE compu- 
tations to quantify the degree of chaotic advection have 
been reported in the literature before, their combination 
for 3D performance quantiflcation of realistic chaotic mi- 
cromixers has to our knowledge not been published be- 
fore. The numerical scheme has the potential to assist 
an experimental optimization since geometrical param- 
eters or fluid properties can easily be changed without 
requiring a new experiment. To demonstrate its applica- 
bility, the scheme is applied to evaluate the parameters 
of the staggered herringbone mixer that lead to improved 
performance. 

In the following sections the lattice Boltzmann method 
which is applied to simulate the fluid flow and the Wolf's 
algorithm from which the flnite time Lyapunov exponents 
are obtained are described. Finally, the numerical results 
and conclusions are given. 



II. SIMULATION METHOD 

The lattice Boltzmann method is used to describe the 
fluid flow. The LB method is a simplified approach to 
solve the Boltzmann equation in discrete space, time and 
with a limited set of discrete velocities.'^S The Boltzmann 
equation, given as 

5J + c-V/ = 17(/), (1) 

represents the evolution of the velocity distribution func- 
tion by molecular transport and binary intermolecular 
colhsions. /(f, c, t) represents the distribution of veloc- 
ities in continuous position and velocity space, r and 
c respectively. In the LB approach the position x at 
which f{x,Ck,t) is defined, is restricted to a discrete set 
of points on a regular lattice with lattice constant Ax, 
implying that space is discretized. The velocity is re- 
stricted to a set of velocities Ck implying that velocity is 
discretized along specific directions. At denotes the dis- 
crete time step. The model we adopt is a D3Q19 model 
which is a 3 dimensional model with 19 different velocity 
directions, k — 0, 1, ISW^ The right hand side of the 
above equation represents the collision operator which is 
simplified to the linear Bhatnagar-Gross-Krook (BGK) 
form.l2£l In a discretized form the BGK operator is writ- 
ten as 

nk = io{r,''{x,t)^Mx,t)). (2) 

Here, uj is the reciprocal of the relaxation time of the 
system controlling the relaxation towards the Maxwell- 
Boltzmann equilibrium distribution f'^{x, t). By consid- 
ering small velocities and constant temperature, a dis- 
cretized third order Taylor expansion of the above equi- 




FIG. 1. A typical example of a SHM geometry as it is used for the simulations. The dimensions of this channel are 
32x64xD/Aa; lattice units, where D depends on the distance between the grooves d and the number of grooves per half cycle 
n. H is the height of the channel and w is the horizontal length of the long arm. We define the height fraction as a = h/ H, 
width fraction as /3 = w/W , and the distance fraction as 7 = d/ D. The wall at a:: = 32 is not shown in the figure, in order to 
view the inside of the channel. 



librium distribution function can be written as 
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where /^'^ denotes the equilibrium distribution function 
corresponding to the velocity vector 4, Cfc the lat- 
tice weights, p is the density, po a reference density, and 
Cg = {l/V3)Ax/At is the speed of sound. is the 

equilibrium velocity of the fluid, which is shifted from 
the mean velocity by an amount g/cu under the influence 
of a constant acceleration g. The evolution of the LB 
process takes place in two steps: the collision step where 
the velocities are redistributed along the directions of the 
lattice and the propagation step by which they are dis- 
placed along these directions. These discrete LB steps 
are implemented by the equation 

fk{x+Atck,t+At)-fk{x,t) = -ujAt{fl\x,t)^fk{x,t)), 

(4) 

which gives the dynamic evolution of the distribution 
function and is referred to as a discretized Boltzmann 
kinetic equation. The macroscopic fluid density is given 
by 



p[x,t) = Po^fk[x,t) 



(5) 



and the macroscopic fluid velocity in the presence of ex- 
ternal forcing is given bj^^H 



u(x, t) 



Po 



p{x, t) 



^fk{x,t)ck~'^g. (6) 



It can be shown by a Chapman-Enskog expansion that 
the macroscopic fields u and p from the above equations 



fulfill the Navier Stokes equation in the low Mach number 
limit and for isothermal systems.l^ In order to simulate 
a fiuid fiow through microchannels, periodic boundary 
conditions are implemented along the z direction (see 
Fig. [I]) and no-slip bounce back boundary conditions are 
imposed at the channel walls. The experimental setup 
of Stroock et al. consists of a number of repeated cy- 
cles (^15) consisting of two asymmetric half cycles each. 
By simulating one full cycle only and applying periodic 
boundary conditions we can reduce the required compu- 
tational effort substantially. If the system is in steady 
state, we do not expect any substantial differences in the 
flow field in different cycles. Therefore, this approach is 
valid. If, however, the experimental system would not be 
periodic, such a simplification would not be allowed and 
the full geometry would have to be simulated. 

We simulate a fiuid which is hydrodynamically sim- 
ilar to water, fiowing inside a SHM with a cross sec- 
tion of 96 /im X 192 /im. The length of the channel is 
1536/^m, but can be varied in order to always accom- 
modate a full cycle of the herringbone structure. For 
computational efficiency we have chosen a lattice res- 
olution of Aa; = 3 /im resulting in a fixed cross sec- 
tion of 32Aa:: x 64Aa; and system lengths of the order 
of 512Aa;. Such a relatively low resolution is sufficient to 
properly resolve the flow field as can be observed from 
Fig. [2] where the velocity field perpendicular to the flow 
is shown at positions within different half cycles (a,b). 
Fig. [2J:) shows the proflle in flow direction. While the 
first two subfigures nicely demonstrate the swirling mo- 
tion of the flow, the third one depicts how the flow pen- 
etrates between the grooves and stays mostly unaffected 
close to the upper boundary of the channel. Previously, 
we have shown that a well resolved velocity field with less 
than 4 percent error as compared to the analytical solu- 
tion can be obtained even for a resolution of 6-8 lattice 
nodes in the case of a 3D rectangular Poiseuille Row^ 
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We further demonstrated that flow over random rough 
surfaces can be well resolved even if the small est obsta- 
cles are only described by 2-4 lattice units In the LB 
method, the kinematic viscosity is related to the discrete 
time step through the expression v — Cs^(l/a; — /S.t/2). 
Since wAi is chosen to be 1 to minimize artefacts due 
to the mid-grid bounce back boundary conditions and 
the simulated fluid has the kinematic viscosity of water, 
V = 10~^m^s^^, this implies for the current choice of 
Aa; that M = 1.5 x 10"^ s and Cg ^ l.lSms"^ When 
the magnitude of g is 1.2 x 10^"^ ms^^ along the z di- 
rection, the average steady state velocity of the system 
is M ~ 6.0 X 10~^ ms~^ which corresponds to a subsonic 
flow. The Reynolds number Re = uLjv of the flow is 
« 1.3, where L = ^/H"^ -\- is the characteristic length 
of the channel. One set of simulations is obtained for g 
being 0.4 x 10^'^ ms^^ which corresponds to Re « 0.4. 

When the flow simulation has reached its steady state, 
P = 1000 massless and non-interacting tracer particles 
are introduced at fluid nodes in the z — Q plane and 
then their velocities are integrated at each time step. To 
calculate the FTLE from particle trajectories a group of 
five particles forms four pairs, with every 5th particle 
placed at the center and the remaining ones being placed 
at the four nearest off-diagonal neighboring LB lattice 
sites. The particle at the center traces a fiducial orbit. 
With this arrangement we are able to follow 800 particle 
pairs by using only 1000 particles. 

The trajectories are obtained by integrating the vector 
equation of motion 



At 



= u{R,), j-l,...,P 



(7) 



where Rj denotes the position vector of an individual 
tracer particle. The velocity u{Rj) is obtained from the 
discrete LB velocity field through a trilinear interpolation 
scheme. 

"Chaotic" systems in general have the important fea- 
ture that two nearby trajectories diverge exponentially 
in time. The rate at which these trajectories diverge can 
be related to the ability of the flow field to create con- 
ditions for chaotic mixing. The Lyapunov exponent is a 
possible measure of the performance of a micromixer as 
it is related to the rate of stretching of the fluid elements. 
It is defined by 



Ano — lim - In — W. 



(8) 



where is the distance between two trajectories at 

time t which evolve from an initial separation 2?(0). Aoo 
gives the value of the Lyapunov exponent as t tends to 
infinity. Due to the finite size of any microfluidic system 
it is not possible to implement this definition in a sim- 
ulation code to study the performance of a micromixer 
Also, when two trajectories separate from each other, 
this definition does not allow to understand the ongoing 



stretching and folding dynamics of the flow. A quanti- 
tative measure of the mixer performance based on the 
Lyapunov exponent can be obtain ed by using the FTLE 
instead of the previous expression.l^SEll The FTLE takes 
the dynamical process more completely into account and 
provides a numerically implementable scheme for quan- 
tification of the performance of mixers. The FTLE is 
defined apl 



_ 1 V{t + St) 

ApTLE-^h,^^, 



(9) 



where t is any particular instant of time and St is a finite 
time after which the FTLE is measured. The same pro- 
cess is repeated over N times, where iV is a large number 
denoting the number of times the FTLE is evaluated from 
trajectories of particle pairs. Amon et al. named the 
FTLE as finite time Lagrangian Lyapunov exponent .'^i' 
The convergence of the average FTLE to the Lyapunov 
exponent for large N is discussed in the paper by Tang 
and Boozer,— 



lim (Aftle)^ 



A. 



(10) 



Wolf et al. suggested a method to calculate the FTLE 
from a set of expe rimen tal data which is well applicable 
to our simulations .'^^'^ The key idea of the method is 
to monitor the distance between particles forming a pair 
and to renormalize it by moving back one of the two par- 
ticles towards the other one in case they have separated 
more than a given threshold distance. The FTLE is then 
computed from the sum of individual Lyapunov expo- 
nents measured between replacements. The method has 
been verified on systems with known Lyapunov spectra 
and exact results have been achieved. Following Wolf's 
approach, we implement the following equation to quan- 
tify the mixer performance on the basis of the average 
FTLE as 



i=0 



n) 



v{t. 



(11) 



where ti is the ith time when a FTLE is evaluated, 
25(^1 + Ti) and 'D{ti) are the distance between particle 
pairs at time step ti + Ti and ti, respectively. Ti is the 
number of time steps which a pair of particles take until 
the next replacement and its magnitude could be differ- 
ent for each measure. N is the total number of replace- 
ments made until time t when (A) n is evaluated. If (A) at 
has a positive and non-zero value the particles separate 
from each other at an exponential rate. These particle 
pairs are initially very close to each other and evolve in 
time. If the separation between the pair is greater than 
a maximum distance, the distance between the particles 
is re-adjusted to the initial distance I?(to). For the im- 
plementation of the scheme, for every particle pair one 
of the trajectories is chosen as the fiducial path, while 
the position of the other particle is replaced if the dis- 
tance becomes larger than the threshold value. This is 
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FIG. 2. The velocity field perpendicular to the flow direction at z = 30 and z — 460 is shown in insets a) and b). The figures 
depict the circulating motion of the flow which is triggered by the asymmetry of the herringbone shaped surface structures. 
Inset c) shows the velocity field in fiow direction at y — 21. Due to the choice of P — 0.66 this value corresponds to a position 
close to the tips of the herringbones in case of the first half cycle and to a case in the center of the long arm in the case of the 
second half cycle. As can be observed the velocity vectors point strongly downwards within the first half of the channel, but 
mostly upwards in the second half. 



schematically represented in Fig. [3] The choice of the 
maximum distance is based on a variation of it together 
with a comparison of the obtained FTLE: if it is chosen 
too large, many tracer particles hit the channel walls and 
do not separate any further resulting in a too low value 
of the FTLE. If it is chosen too small, the tracers do not 
have a sufficient amount of time to separate sufficiently 
so that an exponential increase of the distance cannot be 
detected. The chosen value H/2 has been found to be a 
good compromise between the two extreme cases. In or- 
der to avoid errors due to orientation one of the particles 
of a pair is placed along the line of separation. In case 
a replacement point cannot be found due to a wall node 
present at the location, a nearby fluid node is selected 
as the replacement point. If even such points cannot be 
found since all surrounding nodes are surface nodes, the 
replacement is postponed until a suitable replacement is 
possible. 

The following section presents how FTLE can be uti- 
lized for an optimization strategy for chaotic micromix- 
ers. As an example, the influence of different parameters 
which directly affect the performance of the SHM is eval- 
uated. These are the ratio of the height of the grooves 




FIG. 3. A schematic representation of Wolf's method. 
is the distance between a particle pair at an arbitrary time ti. 
If the distance is greater than a maximum distance then one 
of the particles is replaced near to the other particle along the 
line of separation. 



to the height of the channel a, the ratio of the horizontal 
length of the long arm to the channel width /3, the ra- 
tio of distance between the grooves to the length of the 
channel 7 and the number of grooves per half cycle n. 
The width of the grooves is kept fixed at 24 /im for all 
simulations. Fig. [l] provides a pictorial representation of 
these parameters. 
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III. RESULTS AND DISCUSSION 

The performance of the SHM is studied by varying four 
geometrical parameters. While keeping all other param- 
eters fixed (7 =0.089, a =0.2, n=10), the width fraction 
(/3 = w/W) is varied within the range of 0.5 and 0.82 
and the distance fraction (7 = d/D) from 0.04 to 0.11. 
Then, the number of grooves per half cycle (n) is varied 
from 2 to 10 and the height fraction [a — h/H) from 
0.125 to 0.343. The optimization of the SHM as pre- 
sented here is meant to demonstrate the feasibility of the 
algorithm only since several optimization studies of the 
SHM are already available in the literature. Therefore 
our optimization is reduced to a limited variation of the 
four dimensional parameter space. This surely leads to 
a local optimum of the geometrical parameters, but it 
is not assured that the global optimum has been found. 
However, as shown below, our parameters corespond to 
the optimal ones found by other groups suggesting that 
our local optimum coincides with the global optimum. 

Fig. [4] depicts simulated (A) n {t) for different width 
fractions /3 = 0.66, 0.71, and 0.82. The Reynolds num- 
ber is kept fixed at Re = 1.3 and {X)N{t) is obtained 
from tracing the trajectories of 1000 particles. From 
Fig. [4] it can be observed that in each case (A)jv fluctu- 
ates before finally converging to a particular value after 
^ 6.0 X 10^ time steps. All further simulations are run 
until the FTLE have thoroughly converged. The effect of 
the geometry can be measured by comparing the average 
of the converged FTLE which is denoted by A. The error 
bars in Figs. [5] to [8] are given by the standard deviation 
of the data from the point where it has converged. 
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FIG. 4. The average FTLE for 1,000 particles is shown 
versus the time steps for Re = 1.3. After 6.0 x 10^ time steps, 
(A) AT is converged and the average of this converged value is 
denoted by A. The value of A depends on the width fraction 
P as it is further investigated below. 

When one of the parameters is varied it is taken care 
that the other parameters remain unchanged, because 
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FIG. 5. The variation of the maximum averaged finite time 
Lyapunov exponent A with diff'erent width fraction /3, for two 
difi'erent Reynolds numbers. The data indicates that the max- 
imum A can be obtained for a width fraction of /3 = 2/3. Error 
bars are given by the standard deviation of the mean value of 
(A)^. 



we are interested in the dependence of the individual pa- 
rameters on the performance of the SHM. Fig. [5] shows 
the variance of A and as such the performance of the 
SHM with respect to the parameter /3 for two different 
Reynolds numbers. Re = 0.4 and 1.3. Larger values are 
difficult to obtain due to the limits of the LB method or 
would include a substantial increase in computing time. 
Due to the symmetry of the mixer geometry, only values 
for /3 > 0.5 are shown. The datasets peak at /? = 2/3, 
which implies that the degree of chaotic advection is max- 
imized for this particular value of the width fraction /?. 
According to the units chosen above this corresponds to 
w = 130 /im. The curves for the different Reynolds num- 
bers depict that changing the driving force of the fluid 
does influence the absolute value of A, but has no influ- 
ence on the general shape of the curve. Similar studies 
of the Re dependence for other geometrical parameters 
and various different driving forces confirm this behavior. 
Our findings are consistent with the original experimen- 
tal work of Stroock et aL^i^ as well as 2D numerical op- 
timizations by Stroock and McGraw.122! The latter study 

_ V((C-(C) 



how the so-called heterogeneity factor / 



of 



a dye concentration C varies with the number of cycles in 
the mixer. In both publications it is found that /3 = 2/3 
generates a maximum swirling motion of the fluid. 

For the next set of simulations /3 is fixed at the opti- 
mized value of 2/3 and the distance fraction 7 is varied 
from 0.04 to 0.11. Since it is observed from the pre- 
vious results that the optimal value is independent of 
the Reynolds number, the following simulations are per- 
formed at a constant Re = 1.3. The average value of 
converged FTLE for different distance fractions is shown 
in Fig. [6j It can be observed that after a moderate in- 
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crease of A with 7, the curve has a sharp peak at 7 = 0.07, 
which corresponds to a value oi d = 105 /zm for the cur- 
rent choice of Ax. Afterwards, A decreases in a similar 
fashion as for small 7, but still at higher absolute values. 
A possible explanation is as follows: when the grooves 
are very close to each other, they create "dead spaces", 
i.e. regions in the micro-channel where the fluid gets 
trapped and cannot move freely. With increasing the 
distance between the grooves, the transverse component 
of the velocity increases, hence "chaotic advection" is en- 
hanced resulting a large value for A. On the other hand, 
if the distance between grooves is too large, the mixer 
behaves like a plain channel without any chaotic advec- 
tion component. The maximum in Fig. |6]is then given 
by the optimal interplay of these two effects. 



rotations and come back to its original position before 
entering the next half cycle. An optimized value for n 
therefore depends on the ratio of optimum rotation to 
the frequency at which the distortions at the end of a half 
cycle occur. Similar to the work presented in the current 
paper, Li and Chen performed LB simulations and used 
tracers to follow the flow field. They, however, quantify 
mixing by computing the standard deviation of the local 
tracer concentration and conclude that SHM with n = 5 
or n = 6 perform best. Even though this result is in 
agreement with our finding, the FTLE analysis clearly 
shows that the channel with n = 5 performs better than 
the one with n = 6. 

0.50 I 1 1 1 1 1 1 1 1 1 1 1 
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FIG. 6. The variation of A versus the distance fraction 7 for 
Re = 1.3 and /3 = 2/3. The FTLE rises with the increase of 
7 until it reaches a distinct peak. After that the data shows 
a decrease, indicating that the optimized performance of the 
micromixer is at a groove distance 7 = 0.07. The error bars 
are given by the standard deviation from the mean. 

After having optimized the values for /3 and 7, the 
number of grooves per half-cycle n is varied from 2 to 
10. It can be understood from Fig. [7] that a variation 
of n has the largest impact on the performance of the 
mixer as compared to /3 or 7. For the current setup, by 
variation of n it is possible to change the value of A by 
a factor of 2.3 as compared to 1.2 for /3 and 1.3 for 7. 
Fig. [7] clearly shows that a staggered herringbone mixer 
with n ~ 5 performs best. The existence of a maximum 
performance in dependence on the number of grooves can 
be explained by the fact that interplay between advection 
in flow direction and the swirling motion needs to be 
optimized for a well performing mixer. If the number 
of grooves is too small, the flow fleld is not sufficiently 
rotated when flowing through a half cycle. The change of 
direction of the swirling motion at the beginning of the 
following half cycle does not result in a relevant distortion 
of the trajectories then. If the number of grooves is too 
large, however, the fluid might perform one or more full 
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FIG. 7. The variation of A with the number of grooves per 
half cycle (n) is shown. It can be observed that the SHM 
with n = 5 performs best. Error bars indicate the standard 
deviation from the mean values. 

The final parameter to be varied is the ratio of the half 
depth of the grooves to the height of the channel a. Fig. [8] 
shows the average value of the converged Lyapunov ex- 
ponents for different a between 0.125 and 0.343. After 
a strong increase of A (a), the data shows a maximum at 
a = 0.25. For the units chosen above this corresponds 
to a groove depth of 24 /ini. For larger a the value of 
A decreases again. Our result is similar to the original 
experimental analysis of Stroock et al.'^ Again, an argu- 
ment can be found for the existence of an optimal value 
for the ratio between groove depth and system height: if 
the grooves are too shallow, they are not able to generate 
the swirling motion required for chaotic advection. On 
the contrary, if the grooves are too deep the flow is not 
able to fully penetrate the grooves. Further, the volume 
between grooves and top surface becomes so small that 
the swirling motion cannot develop anymore. 

In this section we have demonstrated that a numerical 
scheme based on the LB method for the flow in com- 
plex mixer geometries together with Wolf's method to 
calculate FTLE from trajectories of passive tracers are 
a powerful tool for the quantiflcation of chaotic mixing. 
Without loosing generality, a limited optimization study 
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FIG. 8. The variation of the average converged FTLE with 
the height fraction (a). The data indicates that the maximum 
FTLE can be obtained for a — 0.25. Error bars indicate the 
standard deviation from the mean values. 



was performed for the particular example of the SHM. 
The optimal parameters a = 0.25, /3 = 2/3, 7 = 0.07 
and n — 5 have been found which is in good agreement 
with known experimental and numerical literature data. 



IV. CONCLUSION 

Mixing at the microscale can be efScient if a large in- 
terface between fluids is provided. This can be obtained 
by passive chaotic micromixers utilizing repeated stretch- 
ing and folding of the fluid interfaces. The performance 
of such mixers depends on the rate at which "chaotic ad- 
vection" of the fluid takes place. For the development of 
efhcient chaotic micromixers it is mandatory to under- 
stand the underlying transport processes as well as their 
dependence on the geometric structure of the microfluidic 
device. In this paper we have demonstrated an efficient 
numerical scheme which allows the quantification of the 
performance of a micromixer. The scheme is based on 
a LB solver to describe the time dependent flow field in 
complex mixer geometries combined with Wolf's method 
to compute FTLE from passive tracer trajectories. We 
have demonstrated the applicability of the quantification 
method by applying it to find an optimal geometrical 
configuration of the SHM, but the scheme should be gen- 
erally applicable to any chaotic mixer. By performing a 
systematic variation of the relevant geometrical parame- 
ters we obtained a set of optimal values which is consis- 
tent with literature data published by others. However, 
those data was obtained from experiments or numerical 
simulations not taking the chaotic nature of the tracer 
trajectories into account. The method presented here, 
however allows a quantification and optimization of the 
mixing performance by investigation of the underlying 
flow profiles. Further work could include the optimiza- 



tion of other chaotic micromixer geometries and the in- 
troduction of multiple fluids to the problem. For the 
latter, the lattice Boltzmann method offers a number of 
possibilities to simulate multicomponent flows and is thus 
a well suited candidate. 
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